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O ' Accurate estimation of licensed channel Primary User's (PU) temporal statistics is important for 
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Dynamic Spectrum Access (DSA) systems. With accurate estimation of the mean duty cycle, u, and the 
mean off- and on-times of PUs, DSA systems can more efficiently assign PU resources to its subscribers, 
thus, increasing channel utilization. This paper presents a mathematical analysis of the accuracy of 
estimating u, as well as the PU mean off- and on-times, where the estimation accuracy is expressed 

^, I in terms of the Cramer-Rao bound on the mean squared estimation error The analysis applies for the 

.. . traffic model assuming exponentially distributed PU off- and on-times, which is a common model in 

K^ ' traffic literature. The estimation accuracy is quantified as a function of the number of samples and 

m 

0> ' observation window length, hence, this work provides guidelines on traffic parameters estimation for 

^\i , both energy-constrained and delay-constrained applications. For estimating u, we derive the mean squared 

CO . estimation error for uniform, non-uniform, and weighted sample stream averaging, as well as maximum 

o 

psj ' likelihood (ML) estimation. The estimation accuracy of the mean PU off- and on-times is studied when 

ML estimation is employed. Besides, the impact of spectrum sensing errors on the estimation accuracy is 
, studied analytically for the averaging estimators, while simulation results are used for the ML estimators. 

X ■ 

^ . Furthermore, we develop algorithms for the blind estimation of the traffic parameters based on the derived 

' ' ' theoretical estimation accuracy expressions. We show that the estimation error for all traffic parameters is 

lower bounded for a fixed observation window length due to the correlation between the traffic samples. 
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On the other hand, the impact of spectrum sensing errors on the estimation error of u can be ehminated 
by increasing the number of traffic samples for a fixed observation window length. Finally, we prove 
that for estimating u under perfect knowledge of either the mean PU off- or on-time, ML estimation can 
yield the same estimation error as weighted sample averaging using only half the observation window 
length. 

I. Introduction 

Spectrum sensing is the cornerstone of Dynamic Spectrum Access (DSA) HI where Secondary Users 
(SUs) search for, and operate on, licensed spectrum that is temporarily vacant. The SUs have to sense 
for the presence of Primary (licensed) Users (PUs) on the targeted spectral bands before utilizing these 
radio resources. The PU channel utilization patterns are stochastic in nature Q. Consequently, acquiring 
knowledge about the PU traffic statistics can improve the performance of SU channel selection algorithms, 
for example ||3], and help in achieving more efficient resource allocation, for example \j4j, in DSA systems. 

A. The Need for Accurate PU Traffic Estimation: an Example 

The multi-channel Medium Access Control (MAC) protocol proposed in fSl is a good example for 
showing the importance of PU traffic parameters estimation. In the proposed MAC protocol, the SUs 
access PU channels opportunistically and sense the presence of PUs periodically. The sensing period 
for each channel is optimized to maximize the expected throughput by minimizing fSl Eq. (1)] which 
quantifies the sensing overhead (denoted by SSOH) and the missed channel access opportunities (denoted 
by UOPP). The optimal sensing period is derived as a function of the PU traffic parameters, specifically, 
the mean PU off-time, 1/Aj, and the mean PU duty cycle, u. We show that when the PU traffic parameters 
estimation enor increases, the performance of the proposed MAC protocol (measured in terms of UOPP 
and SSOH) deteriorates. The results of the investigation on the impact of the PU traffic parameters 
estimation error on this MAC protocol are presented in Fig. [T] We observe that as the deviations between 
the actual and estimated (i) mean PU off-time (Fig. 1 1(a) i and (ii) mean PU duty cycle (Fig. 1 1(b)] ) increase. 



the level of sensing overhead and missed opportunities, SSOH-i-UOPP, increase. For example, even when 
the estimation error in u is only 15%, the resulting SSOH-i-UOPP exceeds the optimal value (i.e. having 
perfect estimates of PU traffic parameters) by almost 10%. Furthermore, we observe that inaccurately 
estimated u has a more profound impact on the performance of the MAC protocol, than inaccurately 
estimated mean off-time. 
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Fig. 1. Impact of inaccurate estimation of PU departure rate for channel i, \j\i, (Fig. |l(a)| i and mean PU duty cycle for 
channel i, Ui, (Fig. |l(b)"| > on the performance of the multi-channel MAC protocol proposed by Kim et al. in (5). The results were 
obtained as follows. With the number of PU channels, TV, set to 2 and sensing time set to 1.6ms, vector A = {A/,i, • • • , A/,jv} 
in Fig. |l(a)| ([/ = {u\, ■ ■ ■ ,un} in Fig. |l(b)) . was shifted from the actual value (by the same factor) positively for the first 
channel and negatively for the second channel. Then, the corresponding optimal sensing period was calculated for the actual and 
erroneous A and U using ||5] Eq. (1)]. The resulting percentage increase in SSOH+UOPP was calculated for the case of having 
erroneous A, (Fig. 1 1(a)") , and erroneous U, (Fig. |l(b)] l. For both figures inter-sampling granularity and maximum inter-sampling 
period were set to 0.5 ms and 0.1 s, respectively. 



B. Related Work 

A large number of algorithms in DSA systems, considering all layers of the communication stack, 
assume perfect knowledge of the PUs' traffic parameters, see for example IH Sec. 2.1], Sec. II- 
B], m Sec. II], d Sec. Ill], QOl Sec. 3.1], [iL Sec. 3.2], HH Sec. II], IHl Sec. II-A], O Sec. III-A]. 
These parameters include the mean PU duty cycle, and the mean PU off-time and on-time. In reality, 
however, DSA systems need to periodically estimate the level of PU traffic before making any decisions 
on PU channel access. As DSA systems often cannot assume any a priori knowledge regarding the PU 
traffic parameters of the accessed channels, blind or semi-blind estimation methods of time-domain PU 
channel occupancy statistics need to be employed. Therefore, the issue of efficient estimation of traffic 
parameters of the PU, considering analytical models of the estimation process, started to gain attention 
from the research community. Recently published |[T5l is a good example of a DSA system where the 
need for the most accurate estimation of PU traffic is essential. Therein, a system which scavenges 
spectrum opportunities in the range of (0,400] milliseconds is introduced and implemented on the TelosB 
mote (TI MSP40 microcontroller, Chipcon CC2420 transceiver) operating on a 2.4 GHz range wireless 
sensor network testbed. One of the components of the designed channel access engine is the channel 



measurement and modeling component. Unfortunately, the paper does not discuss how to design such a 
module. Moreover, two designed channel access strategies: (i) Contiguous Secondary User Transmission 
and (ii) Divided Secondary User Transmission, rely strongly on the knowledge of PU traffic (denoted 
as "whitespace probability density function" by the authors). All PU traffic profiles were artificially 
generated beforehand and known to the channel access engine, which is an unrealistic assumption. 

The most notable results dealing with analytical estimation of PU time-domain traffic parameters 
can be found in ||5l, llT6l - |[T8l . For analytical tractability, all considered works assume that PUs have 
exponentially distributed off- and on-times. In Q maximum likelihood estimation was adopted for 
estimating the mean PU off-time while sample stream averaging was used for estimating the mean 
PU duty cycle. Meanwhile in llT6l . Bayesian estimation was proposed for estimating the mean PU off- 
and on- times. Uniform traffic sampling was assumed for both ||5l and |[T6l . On the other hand, the 
authors in IITtI . lITSl . using the notion of Fisher information, derived optimal traffic sampling schemes 
for estimating the mean PU off-time. They argued that for a fixed channel observation window and a 
fixed number of samples, random sampling outperforms uniform sampling. However, perfect knowledge 
of the mean PU duty cycle was assumed. Besides, no closed form expressions for the accuracy of the 
estimated mean PU off-time was derived and different random sampling schemes were evaluated only via 
simulations. Unfortunately, in all aforementioned works O, |[T6l - |[T8l . the estimation accuracy, measured 
in terms of the mean squared error (MSE) in the estimated parameters, was not quantified in a closed form. 
Moreover, the impact of spectrum sensing errors on the estimation accuracy was not studied analytically. 
In |[T9l the authors derived the bounds on the accuracy of the joint estimation of the arrival and departure 
rates of PUs. However, the authors assumed that the PU traffic is observed continuously, which is an 
assumption that is far from being practical as the PU traffic is sampled according to a discrete sampling 
process. Also, just like in earlier works, the impact of spectrum sensing errors in ||T9l was not considered. 

In the context of our work we need to refer to other studies on PU traffic estimation. Specifically, ll20l 
followed a different approach for estimating the PU channel usage statistics (i.e. its complete distribution) 
by using a combination of statistical distance metrics: kernel density estimation, goodness-of-fit testing 
(utilizing the Kolmogorov-Smimov test), and KuUback-Leibler distance. To increase the complexity of 
the problem, cooperation between spatially separated DSA nodes was considered resulting in node-to- 
node variances in PU traffic observations. Unfortunately, no closed-form expressions for the PU traffic 
distribution estimation accuracy were presented. Only a heuristic estimator (in the form of the algorithm 
presented in Table I of ||20l ) was used. The proposed heuristic estimator is based on an example of a 
utility function. Moreover, the impact of spectrum sensing errors was not considered (however, errors 



due to fading and propagation characteristics were included). 

Finally, |[211 . Il22l considered the estimation of the PU channel state through randomized channel 
probing. These papers modeled the PU state estimation problem as an exploration/exploitation problem 
and based the analysis on multi-armed bandit formulation. The difference between these two papers lies 
in system model assumptions and new features that have not been considered in earlier works on multi- 
armed bandit problems for DSA, i.e. lIlTI considered spectrum sensing errors, while ||22| considered PU 
state/channel fading correlation. We need to emphasize however that PU state estimation in ||2TI . ll22l 
has the following limiting features: (i) PU channel state estimation reduces to one parameter only (on 
or off time), (ii) the estimation process requires network feedback, e.g. via ACK/NACK, and (iii) the 
estimator does not collect statistics on the PU channel usage. 

C. Our Contribution 

In this work, we first consider the problem of estimating the mean PU duty cycle, u. We derive the 
estimation MSyJ in u when sample stream averaging with uniform sampling is used. We extend our 
work to include non-uniform sampling as well as weighted averaging with uniform sampling. Moreover, 
we propose estimating u using maximum likelihood estimation under uniform sampling, and derive the 
corresponding Cramer-Rao (CR) bound which provides a lower bound on the estimation error for unbiased 
estimators employing uniform sampling. Regarding the mean PU off- and on-times, we derive the CR 
estimation error bounds for both parameters under uniform sampling, and present the corresponding 
maximum likelihood estimators. All of the estimation error expressions presented in this work are 
formulated as functions of the total number of samples, which serves as a guideline for energy-constrained 
applications where the energy budget for sampling, and hence the total number of samples, is limited. 
We also quantify the relationship between the estimation error and the length of the observation window. 
This is important for delay-constrained applications, and when non-stationary traffic is considered, as it 
shows the compromise between the delay in learning the PU traffic parameters and the estimation error 
in the parameters. Besides, the effect of spectrum sensing errors on the estimation accuracy is studied 
analytically for the averaging estimators, while simulation results are used for the ML estimators of 
u, and the mean PU off- and on-times. Finally, we use the resulting expressions to design algorithms 

'in parameter estimation literature, the MSE is often used as a metric for tfie estimation accuracy for a number of reasons. 
The MSE is an intuitive metric that describes the average squared deviation of the estimated parameter from the actual value of 
the parameters. Moreover, the MSE accounts in the same manner for both positive and negative deviations. Finally, the MSE 
metric is mathematically tractable and can often be expressed in closed form, as opposed to the mean absolute error, or the error 
probability. Closed form expressions have the advantages of providing intuition regarding the results, and enabling incorporating 
other mathematical tools to analyze the results. 



for the blind estimation of the PU traffic parameters under a variety of constraints, and compare their 
performance against the derived theoretical bounds. 

The paper is organized as follows. Section |II] presents the system model considered in this work. 
Expressions for the MSE in estimating the PU traffic parameters are derived in Section Ull] (duty cycle) and 
Section |IV] (off- and on-time rate parameters). Two practical algorithms for estimating the PU duty cycle 
and mean off- and on-time are presented in Section |Vl while numerical results are given in Section |Vl] 
Finally, Section IVIII concludes the paper. 

II. System Model 

Following the model introduced in Q, we consider a single channel that is licensed to a single PL|3. 
The PU traffic is assumed to be stationary over a sufficiently large time window with exponentially 
distributed off- and on-timefl The probability density function of an exponentially distributed random 
variable, x, is given as |[24l Eq. (3.15)] f\{x) = \e~^^ , for x > and f\{x) = 0, otherwise, where 
A is denoted by the rate parameter. With A = A/ and A = A„, f\{x) denotes the distribution of PU 
off- and on-times, respectiveljQ. The mean PU off- and on-times are equal to the reciprocal of Aj and 
A„, respectively. Besides, the duty cycle u of the PU can be calculated as |[24l Sec. 11.3] u = ^ ' . 
Hence, A/, A„ and u are inter-dependent, where estimating any two of the three parameters is sufficient 
to completely estimate the PU traffic parameters. 

In order to estimate the traffic parameters, the PU channel is sampled in order to acquire data regarding 
the state of the PU (on or off). For the considered system model, denote the total number of samples by 
N. Denote the PU traffic samples by the vector z = [zi,Z2,- • • , z^] where z„ is the nth traffic sample, 
and Zn = 1 if the PU is active and Zn = 0, otherwise. Moreover, in the proposed model, we consider the 
general case where the spectrum sensing process is prone to errors. The sensing error is modeled in the 
form of false alarm and mis-detection probabilities, denoted by Pf and Pm, respectively. The sensing 
error is assumed to be independent for different traffic samples. The estimated PU traffic samples are 

^Note that this, and other assumptions of Q, like, e.g., (i) introduction of (collaborative) spectrum sensing, (ii) listen-before- 
talk policy, (iii) scheduling of quiet periods, (iv) availability of the control channel, are standard and axiomatic in the DSA 
literature. Therefore our results are a natural extension of a well-established path in DSA research. 

^We are considering a continuous model as it is more general than a discrete one, encapsulating the discrete traffic case. 
Furthermore, discrete PU traffic models impose an implicit synchrony between the PU and DSA networks. This requires a priori 
knowledge of the PU properties, e.g. guard intervals or pilot symbols. An example of such operation is in 1231 . where the DSA 
system operates following the slot boundaries of a GSM system. To avoid such constricting requirements, we made as little 
assumptions on the PU properties as possible. 

'*Note that the assumption on the exponential distribution of off- and on-times is common in DSA literature, e.g., see recent 
examples of I25I - I27I : see also recent papers confirming the exponential distribution of time-domain utilization of certain 
hcensed channels JU, ||28l, (29) • 
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denoted by the vector z = [zi,Z2,- " , ^Af] where z„ is the 77,th estimated traffic sample. It follows that 
Zn = 1 if Zn = 1 and no mis-detection error occurred, or Zn = and a false alarm error occurred. 
Similarly, z„ = if 2;„ = 1 and a mis-detection en^or occurred, or z^ = and no false alarm error 
occurred. Furthermore, the inter-sample times are given by the vector T = [Ti, T2, • • • , T^-i] where T„ 
denotes the time between samples Zn and Zn+i- Finally, the total observation window length is denoted 
by T, where T = Zn=i Tn- 

Denote the PU state transition probability by Vvxyit), which corresponds to the probability that the 
PU state changes from state x to state y within time t, where {x,y} = denotes that the PU is idle 
while {rr, y} = 1 denotes that the PU is active. The PU state transition probabilities were derived in ||5] 
Sec. 6.1] as 

1 — u + ue " , x = 0,y = 0, (la) 

l-Proo(i), x = 0,2/ = l, (lb) 

-A if 

u + {l — u)e'~^, x = l,y = l, (Ic) 

l-Prn(i) x = l,y = 0. (Id) 

In this work PTxy{t) is later used to derive the MSB in the estimates of u, Aj, and A„. 

As remarked in Section H estimators of u and Aj are analytically described in closed form in ||5], 
IIT6I - IIT8I . However, a measure of the estimation error in u and Xj, was not given, noting that in ^ 
Sec. 6.2] only the asymptotic confidence interval for the estimates of u and A/ was presented. In the 
following sections, we propose new methods to estimate u and we derive the MSB in the estimates of 
u, A/, and A„. 

III. Estimation of the Primary User Duty Cycle u 

In this section, we analyze different methods for estimating the duty cycle, u, of the PU. We first 
present an estimator based on averaging the traffic samples, labeled the averaging estimator, similar to 
the estimator presented in ||5], llT6l - l|T8l . In addition, we modify the estimator to the general case where the 
PU traffic samples are not uniformly sampled. Furthermore, as we observe that the estimation accuracy 
is bounded by the sample correlation, we propose two different estimation methods to alleviate the 
correlation effect. The first method is based on the weighted averaging of the traffic samples, labeled the 
weighted averaging estimator, and the second method is based on maximum likelihood (ML) estimation. 
For all three estimation methods, we derive expressions for the MSB in the estimates. Moreover, we 
derive the CR bound on the estimation error when using uniformly sampled traffic samples. The MSB 



expressions are presented as functions of the number of samples and the observation window length to 
serve as guidelines for traffic estimation in energy-constrained and delay-constrained systems, respectively. 

A. The Averaging Estimator under Perfect Sensing 

The averaging duty cycle estimator, Ua, is defined as IH Sec. 6.1] 

1 ^ 

"^0. = J^^Zn- (2) 

n=l 

We first consider the case where the spectrum sensing errors can be ignored, i.e., Pf = P^ = 0, hence, 
Zn = Zn'^n. The impact of spectrum sensing errors on the estimation en^or is presented in the next section. 
1) The MSE in Ua-' The MSE in Ua for N samples can be defined as 

Vu^,N = E[ul]-u', (3) 

where E[-] denotes the expectation. The intuition behind ([3]) is as follows; the expectation is calculated 
over all possible values of Ua resulting from all 2^ permutations of the traffic samples vector z. Define 
Z as a. vector containing all 2^ permutations of z with Z„, n G {1,2,- •• ,2^}, defined as the nth 
element of Z. Furthermore, define Zn^m, m G {1, 2, • • • , N}, as the mth traffic sample of Zn, and define 
(n = ^m=i ^n,m, 1-6-. the Summation of all traffic samples of Zn. Then, substituting ([21) in Q yields 

2 

iV2' 



Vu.,N = ^ECnPr(^ = Zn\T) - u\ (4) 



n=l 

where Pr(2 = Zn\T) denotes the probability of observing PU traffic sample sequence Z^, for a given 
vector of sampling times T. We then have the following corollary. 
Corollary 1: The MSE of Ua is given as 

/I \ r, /-, N N-lN-ii+j-l 

Uil — U) luil — U) yr^ JT-^ -i-r '^fc^/ 

Proof: See the proof of ([19]). D 

From Corollary 1 we obtain the subsequent corollary. 



Corollary 2: The decrease in the MSE in Ua with each extra sample, D^^^at+i, is given as 

_ m+1) u{i-u) 

(A^ + l)2 "°'^ (iV+l)2 

(N-l N \ 

i+^E n ^"^ ■ (« 

i=0 k=N-i } 

Proof: Via elementary algebra. D 

Corollary |2] as we will show in Section |V] proves important in designing adaptive algorithms for the 
blind estimation of u. 

Remarks: The rightmost term of (|5]l represents the increase in the estimation error caused by the 
correlation between the traffic samples. As T^ tends to infinity, this term tends to zero, hence, Vu^,N 
approaches "^ ^^' , which is the MSE in estimating the duty cycle of an uncorrected traffic sample 
sequencqj. This is attributed to the fact that the inter-sample time becomes large compared to the mean 
off- and on-times of the PU, hence, the coixelation between the samples vanishes. The estimation error 
is a function of the traffic parameters, the number of samples, and the inter-sample time sequence. The 
optimal inter-sample time sequence that minimizes the estimation error for a given number of samples 
and a fixed total observation window length is derived in the next section. 

2) The Optimal Inter-Sample Time Sequence for Minimizing the MSE in Ua: In this section, the MSE 
in Ua is shown to be convex with respect to the inter-sample time sequence, T. The optimal T, denoted 
by T* , is derived, and the corresponding expression for the MSE in Ua is presented. Expression dSjl is 
proven to be convex by showing that the Hessian of Vu^^n {T), denoted by V^T^^^at (T), is positive- 
semidefinite |[3T1l . The proof of convexity is given in Appendix iBl 

'Note that "^ ~^' is the variance of a binomial distribution normalized by N'^ where the probability of success is set to 
u Oni Ch. 4]. 
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The problem of minimizing Vu^^n (T) with respect to T can be written as: 

u(l 



minimize T4„,7v(7') 



u) 2u{l — u) 

N~lN-ii+j~l 

j=i j=i k=j 



subject to - r„ < 0, 71 = 1, 2, • • • , A^ - 1; 

N~l 



(V) 

(8) 
(9) 



n=l 



The optimization problem can be solved by Lagrangian duality BTl where the Lagrangian function can 
be expressed as 

AT-l 



Lv (T, v,ij) = Vu,,N (T) - ^ VkTk 

fc=i 



(10) 



.fc=i 



where v = [vi,V2,- " ^ ^^Af-i] is the vector of the Lagrangian multipliers associated with inequalities dSjl 
and /i is the Lagrangian multiplier associated with (|9]l- As the optimization problem is convex, and the 
objective and constraint functions are differentiable, the optimal inter-sample time sequence, T*, satisfies 
the Karush-Kuhn-Tucker (KKT) conditions: 

(11a) 



T < 

-1 

J2Ti:-T = 0, 
= 0,<>0, 

tVK„,7v(r„*) -< + /.* = 0, 



-n 
N-1 



k=l 

v*T* 



(lib) 
(lie) 
(lid) 



where nG{l,2---,iV — 1} and the superscript * signifies optimality. Expression (llldl i can be expressed 



as 



-2u{l~u)\f ^— -vn 



mu 



Z-^i=l L^j=n i.lk=i 



-A/T,!/« 



f* + /x* = 0. The solution of the convex problem is first 



presented for the special case where T^ > 0, t;* = Vn, i.e. when all samples have non-zero spacing in 
time. The solution is later expanded to include cases where T* = Q for a set of n, i.e. samples coincide 
in time implying that the sample is weighted by the number of coinciding samples. 



For the case where T* > {) Vn, T* is given as follows. Denote F* = e" 



and r* 



", then 



11 



n 



T* for n 



1 and n = N 



1, and T* = T^, otherwise, where 

r* 



1 



2r; + (iv - 3) Tfc* = r, 



u 



T> (Af-3)-^log2. 



(12a) 
(12b) 
(12c) 



Equation ( |12a| ) is derived by simultaneously solving (llldl ) for n = 1 and n = 2. Equation ( |12b| ) is derived 
from condition (lllbl l. Condition (I12cl i is derived by setting T^^ = in (I12al ) and (I12bl ) and solving for 
T. Expressions (I12al ) and (I12bl ) can be shown to satisfy (lllal i- dlldl ) but the proof is omitted for brevity. 
The solution for the optimization problem implies that as the length of the total observation window, 
T, increases, the optimal inter-sample time sequence approaches uniform sampling. As T decreases, T* 
remains uniform for samples 2 to A^ — 1. However, the first and last inter-sample times are equal in 



length, and shorter in length than the rest of the inter-sample times. If T is decreased to 



{N-3)u 



log 2, 



Tj* and T^_i approach zero, i.e. the first two samples and last two samples coincide. This implies that 
the number of samples is decreased to A^ — 2 and the first and last samples are weighted by two. 



For the case where T < 



{N~3)u 
Af 



log 2, T* can be equal to zero for n G /C where /C = {1, 2, • • • ,k — 
2,k-l,N-k + l,N -k + 2,--- , A^ - 1} and 1< A; < [f J . T* can be found by solving (ITTa-dTTdl) 
where T* = for n € JC and T* > 0, < = 0, otherwise. T* is derived as T^ = for n e )C, T* = T* 
ior n = k and n = N — k and T* = T^, otherwise, where 

r* 



° ~ A; (1 - r*) ' 

2T* + NT^ = T, 

- u , /c + 1 ^ ,% It , 
A^— log —— <T< N— log - 
Af k Af k 



1' 



(13a) 
(13b) 
(13c) 



where N = N — 2k — 1. Equation (I13al ) is derived by simultaneously solving (llldl ) for n = /c and 
n = k + 1. Equation (I13bb is derived from condition (lllbl l. The lower bound in (I13cl ) is derived by 



setting T* = in (I13al ) and (I13bb and solving for T. The upper bound in (I13cl ) is based on the condition 
'^k-i — where the expression for u* can be derived for n G /C from (llldl ). using (I13al ). as v^ = 



2u{l-u)Xf 



{k — n) 



i-r; 



n 



. Again, (I13al ) and (I13bb satisfy (lUaMUdl ) Vn, and the proof is omitted 
for brevity. The solution for the optimization problem implies that if T falls in the boundary expressed in 
(I13cl ). then the first and last k— 1 samples are omitted, and the kth and (A^ — k)th samples are weighted 
by k. Again, the middle N — 2k — 1 inter-sample times are uniformly sampled and T^ = T^_j^. 
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For the special case where < T < {N — 2k + l)^log ^t^ and k = [^J, all inter-sample times 
decrease to zero except for the middle interval, for A^ even, or the middle two intervals, for N odd. 
Accordingly, for even N, T* = ioi n y^ ^ and T^^^ = T, whereas for odd N, T* = for n ^ 
{^, ^} and T %_, = T^i = T/2. This can be shown to satisfy the KKT conditions fTTah- fTTdl) . 

2 2 

Using the optimal inter-sample time sequence, the lower bound on the estimation error for the averaging 
estimator, Ua, can be derived by substituting (I13al ) and (I13bl ) in (JS) yielding 

_ 2u(l - u) 



^ 'N^ri + k{k-i){i-r, 



*\2 



2 (i-r*)2 

r*(iv-2fc)(i-r*: 
+ (i-r*)2 '' ^'^' 



for < A; < [-yj where k is chosen to satisfy ( |13c|) and T^ is found by solving ( |13ab and ( |13br 



Remarks: The optimal inter-sample time sequence, T*, is non-uniform where the first and last inter- 
sample times are shorter than the rest of the inter-sample times. Moreover, generally, the first and last 
samples have higher weights compared to the rest of the samples. This is attributed to the fact that the 
first and last samples are at the edges of the traffic sample sequence, hence, they have lower correlation 
with the rest of the traffic samples. However, as the total observation window length increases, the impact 
of the first and last samples on the overall estimation error decreases, and the gap between the estimation 
error for the optimal non-uniform sampling sequence and the uniform sampling sequence diminishes. 
Furthermore, T* is a function of u, the very same parameter that is to be estimated, as well as a function 
of the mean PU departure rate, Xf (or the mean PU arrival rate, A„, as Aj equals nA„/(l — u)), which is 
not necessarily known by the traffic estimator. Hence, T* cannot be known a priori, yet T* can be used 
as a guideline in algorithm design for the blind estimation of the traffic parameters. For instance, apart 
from 'weighting' the first and last samples, T* is found to be an almost uniformly sampled sequence. 
Besides, the error expression given in ([14) serves as a lower bound on the MSE in estimating u using 
averaging for any inter-sample time sequence. 

3) The Averaging Estimator under Uniform Sampling: The work in ||5], |[T6l - |[T8l considered estimat- 
ing u by averaging uniformly sampled traffic observations. In this section, we derive the MSE in Ua under 
uniform sampling, denoted by Vu^^^n- With constant inter-sample times, T„ = -j^^ = T^, VT„ G T. 



•^For k = [fj, for even N, 14*„,^r = ^"<^^") UL + [N^{e-^i'^/'' + 1) - 2JVJ /4j, and for odd N, V^^ 
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Substituting in ^, the MSE can be written as 



2n(l-u)/iV ^' i,,, .A 



_ 2n(l - u)ru{r^ - NjTu - 1) - 1) 

iv2(i-rj2 

where T^ = e ^. Again, the leftmost part in the second equation of dTSl ) accounts for the increase in 
estimation error caused by sample correlation. Intuitively, when the sample correlation is high, increasing 
A^ in a fixed time window leads to an insignificant change in the estimation error. Formally, we obtain 
the following corollary. 

Corollary 3: For a fixed observation window length, as the number of samples increases, the MSE 
error in estimating u for uniform sampling approaches an asymptote Vu^^^l, where 

14„„,i = lim Vu^^,N = ^^^iipl {e-^ + r? - 1) , (16) 

where ?? = — ^. 

Proof: Via elementary algebra. D 

Note that Vu^^^^l tends to as the observation window length is increased. Using Corollary |3] the 
number of samples, N, can be chosen such that the resulting error is above the asymptotic error (IT&i by 
a factor /3. Then N can be evaluated by solving Va^^^N = /314„„,L- 

Remarks: When estimating u by averaging uniformly sampled traffic samples, the estimation error 
is lower bounded. The lower bound is caused by sample correlation and can only be eliminated by 
increasing the total observation window length. 

B. The Averaging Estimator under Imperfect Sensing 

The analysis presented in Section IIII-AI is extended here to include the effect of spectrum sensing 
errors on the estimation error. Introducing spectrum sensing errors to the averaging estimator expressed 
in Q causes the estimator to become biased. The expected value of the estimator can be calculated as 
E [ua] = Pf {1 — u) + u {1 — Pm), where the expectation is calculated over all possible values of Ua 
resulting from all 2^ permutations of the estimated PU traffic samples vector z. Thus, the duty cycle 
can be calculated from E [ua] where u = {E [ua] — Pf)/{^ — Pf — Pm)- Accordingly, we define the 
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following unbiased estimato 

1 



(17) 



l-Pf-Prr 



n=l 



yPfjPm : Pf + Pm / 1- Define Z as a vector containing all 2^ pemiutations of z with Z„, n G 
{1, 2, • • • , 2^}, defined as the nth element of Z. Furthemiore, define Zn,rn, ?ti G {1, 2, • • • , N}, as the 
mth traffic sample of Zn- Thus, the MSB in Ua,s for N samples, denoted by ^s^^.at can be expressed as 

K„..,N = 5^52pr(i = i„|r)-t/2, (18) 



where Sn = i^p^^p. 



n=l 

N 



Theorem 1: The MSB in n^ ^ is given as 



. We then have the following theorem. 



4 = 1 j = l fc=J 

^P^.(l-P^.) + (l-n)P/(l-P/) 

Proof: See Appendix lAl D 

Remarks: Comparing (|T9] l with dSjl, it is clear that the rightmost term of the right hand side of 
( fT9l ) models the increase in the estimation error caused by the spectrum sensing errors. Moreover, the 
leftmost term of the right hand side of (|T9] l accounts for the estimation error caused by the sample 
correlation. Furthermore, unlike the impact of sample correlation on the estimation error, the effect of 
the spectrum sensing errors on the estimation error can be asymptotically eliminated by increasing A^. 
Besides, the increase in the estimation error attributed to the spectrum sensing errors is not a function 
of the inter-sample time sequence, T. Accordingly, Vu^ ^,n is convex with respect to T, and the optimal 
T that minimizes the MSB in Ua^s is the same as T that minimizes the MSB in Ua that is derived in 
Section IIII-A2I 

^Note that the estimator Ua,s is not defined for the special case of Pf + Pm ~ 1. For Pf + Pm = 1, the denominator of the 
proposed unbiased estimator equals zero, and the expectation of the biased estimator can be expressed as E [ua\ = Pf, which 
is independent of u. Hence, both estimators fail to estimate u. On the other hand, note that Pf + Pm > 1 does not correspond 
to any relevant practical sensing method. Typical values for the probability of false alarm and mis-detection are Pf ^ 0.1 and 
Pm ^ 0.1, respectively, e.g., ^ Sec. VII-C], ||33l Sec. VI-A], (H, 135], ||36l Sec. 6.6], [HI Sec. IV-A]. 
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C. The Weighted Averaging Estimator under Perfect Sensing 

In the previous sections, the PU duty cycle, n, is estimated using equal weight averaging of the 
channel samples. The optimal inter-sample times were found to reach zero for some samples implying 
that weighting might improve the estimation accuracy. In this section, we propose a new estimator that 
averages weighted traffic samples to decrease the estimation error by alleviating the effect of sample 
correlation. For analytical tractability, uniform sampling is assumed with a constant inter-sample time 
denoted by Tc. 

We first present the special case where the spectrum sensing errors can be neglected, that is, Pf = 
Pm = 0, and Zn = Zj^n. The effect of spectrum sensing errors on the estimation error is investigated in 

N 

the next section. The estimator is defined as u^ = ^ WiZi, where Wi is the weight of sample Zj. Then 



i=l 



N 



E[uw\ = u"}^ Wi, thus, for the estimator to be unbiased, that is i?[n^] = u, the weights must satisfy the 

i=l 

N 

condition ^ t(;j = 1. 

i=l 

1) The MSE in n^,.- The MSB in n^, can be written as 



T4„,iv = E [{u^ - E[uw]y 



E 



E 



N 




- N 


V" 


)_^ w,z, - 


-E 


y^wiZi 




i=l 




.i=l 


/ 






N 



u I 2_^Wi 



N 



y^ tff + 2 ^ WiWjE[ziZj\ 

i<j 



i=l 



/ N 

\i=i 



(20) 
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The expression E[ziZj] in (l20l) represents the correlation between Zi and Zj denoted by Rij. Consider 
RiA+j, Vj > 1, then 

Ri,i+j = E[ziZi^j] = Prjzj = l,2i+j = 1} 
= Pr{zi = l,Zi+j^i = l}Prii(Tc) 

+ Pr{z, = l,z,+j_i=0}Proi(r,) 

= Ri,i+j-~l Pl'll(T'c) 

+ (n-i2i,,+j_i)Proi(r,). (21) 

The initial condition for the recursive equation (|2Ti is i?j j = -^[ziZi] = u. Thus, solving equation (|2TI) 
yields RiA+j = uFc + ii^(l — Pc), where Pc = e ^^ . Since the traffic samples have a constant mean and 
the correlation function is only related to the time difference between the samples, the samples follow a 
wide-sense stationary processlj and Ri^i+j = R[j],'ij > 0. Substituting Rij in (l20l ) yields 

N N-1 N-j 

i=l j=l i=l 

N N-1 N-j 

i=l j=l i=l 

xu(l-u). (22) 

Note that (l22l ) matches (fTSl) if constant weighting is assumed, i.e., Wi = l/N, i G {1, 2, • • • , N}. 

2) The Optimal Weighting Sequence: The optimal weighting sequence that minimizes the MSE in Uw, 
denoted by w* = [w*, Wg, • • • , w*j^]^ , is derived in this section. According to the orthogonality principle, 
v^* satisfies E[{u — Uu,)zj] = 0,Vj G {1,2, ••• ,N}. Hence, v^* = R^^q„^, where R = [vij]Ny.N, 
Tij = R[\i — j\], and q^^ is the cross correlation vector ii^[uz]. Since the traffic samples follow a wide- 
sense stationary process, R is both symmetric and Toeplitz. The cross correlation vector, q^^, can be 
written as q^^ = E[uz] = u{E[zi], E[z2], ■ ■ ■ ,E[z]\[])'^ = [u'^,u'^,--- ,ti^]^. Normalizing the weights 
to get an unbiased estimator, the optimal weighting sequence is given by w* = ^t^i' , where c is 

*It is stated in (5) that the traffic samples follow a semi-Markov process. But given the condition of using a constant inter- 
sample time, it turns out to be a wide-sense stationary process. 
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a constant vector [1, 1, • • • ,1]^ G M^^^. Hence, the optimal weighting sequence can be expressed as 

The MSE in u^ when the optimal weighting sequence is used can be found by substituting w* in 
(l22l) yielding V^ j^ = ^L"^^\ ^p\ Due to the correlation between the samples, for a given observation 
window length, the MSE in u^ approaches an asymptote, denoted by V-* ^, as N increases. F-* ^ can 
be calculated as follows 

Vl,L = lim Vl^^ = ^^^1^. (23) 

Remarks: The derived optimal weighting sequence dictates that the first and last samples have to be 
multiplied by higher weights than the rest of the samples. This is because they are less correlated to the 
rest of the traffic samples and hence, hold more information. The optimal weighting sequence, however, 
depends on the actual value of u which is obviously not known a priori. Thus, V^ ^ serves as a lower 
bound on the estimation error when weighted sampling is used. Besides, the accuracy of the estimate 
can be improved in an iterative manner where u^ can be used to calculate an estimate of the optimal 
weighting sample, then the resulting weighting sequence can be used to improve the estimate of u. 

D. The Weighted Averaging Estimator under Imperfect Sensing 

In this section, we consider the performance of the weighted averaging estimator considering spectrum 
sensing errors. The spectrum sensing errors cause the estimator to become biased, akin to the averaging 
estimator case presented in Section ITlI-B I Accordingly, the bias can be eliminated by defining the following 
estimator 



"-'^ - l-Pj-Pm 



N 

-Pf + y^^wjZj 



i=l 



(24) 



where wi is the weight of the estimated traffic sample Zj. Hence, E[uw.s\ = u under the condition 
Z]i=i Wi = \.To quantify the 
traffic samples, Ri^i+j, where 



Si=i "^i = 1- To quantify the MSE in u-^i^s, we start by evaluating the correlation between the estimated 



= Pr{zj = l,Zi+j = 1} 

= P]PT{Zi = 0, Zi+, = 0} + P/(l - P„) 

X Fr{zi = 0, Zi+j = 1} + Pr{zi = 1, Zi+j = 0} 

+ (1 - PmfPr{z^ = 1, Zi+j = 1}. (25) 
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Recall from Section BlI-C II that Pi{zi = l,Zi^j = 1} = Ri^i+j. It follows that Prjzj = l,Zj_|_j = 0} = 
u - Ri,i+j, Pr{zi = 0,Zi+j = 1} = u- Ri,i+j, and Prjzj = 0,Zi+j = 0} = 1 - 2u + Ri,i+j. Hence, 
substituting in ( [25] ) yields 



R^ 



i,i+j 



^i,i+j\^ ^rr 



Pf 



+ 2uPf{l-Pra-Pf) + Pf. 



(26) 



Since, as in Section IIII-C1[ the estimated samples follow a wide-sense stationary process, then Rt^t+j 
R[j],^j > 0. Hence, the variance of the proposed estimator can be expressed as 



Vf 



Uu, s,N 




X E 



+p 



{1-Pf- Pm) 



N-lN-j 

+2 I] 5] ^^^i+J^iJ] - '^PfE[H + P] 
i=i i=i 



(27) 



Substituting ^ in (|271) and replacing E[zi] and E[zj] by {I - Pf - P^)u + Pf yields 

'uPrn{l - Pm) + (1 " u)Pf{l - Pf) 



K„ 



.N 



{l-Pf- P™)2 



7 

N N-l N-j 

i=l j=l i=l 




(28) 



In order to derive the optimal weighting sequence that minimizes the MSE in Uuj,s, we apply the 
orthogonality principle, i.e., E[{u — Uw,s)zj] = 0, Vj € {1, 2, • • • , N}. By solving the normal equations 
Rw = qui, where R G M^^^ is the autocorrelation matrix pf^s;^ 



vector [wi,W2r'' ,wn] , and quz is the cross correlation vector E 



^_p _p E[zz ], w is the weights 

Pf 



u + 



the normalized optimal weights as wl 



wl 



-, and Wj 



l-Pf-Pr,. 

1-r, 



, we derive 



This is the same result as for the weighted estimator assuming perfect spectrum sensing. Substituting the 
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optimal weights in (1281 ) yields 



u{l-u){l + T 



C 



N{i - r,) + 2r, 

uPm{l - Pm) + (1 - u)Pf{l - Pf) 



+ 



(l-Pf-Pr 

2 + (iV-2)(l-r ^2 



'^ (29) 



(Ar(l _ r,) + 2re)2 ■ 

Note that the rightmost part of the right hand side of (|29] represents the increase in the estimation error 
attributed to the spectrum sensing errors, while the leftmost part represents the estimation error under 
perfect sensing as derived in Section IIII-C2I It follows that, for a given observation window length, the 
MSE in n^ s approaches an asymptote, denoted by V^ i, as N increases, due to sample correlation, 
where V^ ^ = limA^^oo Vi ^ = ^^%^. 

Remarks: The optimal weighting sequence is not affected by the spectrum sensing errors, as the 
sensing error is assumed to be independent for the different traffic samples. Moreover, as A^ is increased, 
the impact of spectrum sensing errors is eliminated where V^ l ~ ^i L- ^^ ^^^ other hand, the 
estimation error caused by the sample correlation serves as a lower bound that can only be decreased 
by increasing the total observation window length. Furthermore, as in the case of weighted averaging 
assuming perfect sensing, the optimal weighting sequence depends on the actual value of u which is 
obviously not known a priori. Thus, V-*^ ^ ^ serves as a lower bound on the estimation error. 

E. ML Estimation of u and the CR bound on the estimation error under Perfect Sensing 

In the previous sections, it was shown that the accuracy of the estimators of u that are based on sample 
stream averaging is limited by sample correlation. Hence, in this section, we propose a more accurate 
estimator of u based on ML estimation. However, the improved accuracy of ML estimators comes at the 
expense of an increase in the computational complexity. ML estimation is used to estimate parameters 
of a statistical model by finding the parameters' values that maximize the probability of the observed 
samples |[38l . The mean squared estimation error of ML estimators is often quantified analytically using 
the CR bound. The CR bound quantifies the minimum mean squared estimation error that can be achieved 
by any unbiased estimator. ML estimators achieve the CR lower bound as the sample size tends to infinity 
when certain conditions are satisfied ||39] Ch. 12]. Accordingly, we present the likelihood function for the 
estimation of u, as well as the corresponding CR bound on the estimation error. The expressions are first 
presented for the special case when the sensing procedure is assumed to be perfect. Then, the likelihood 
function is modified to account for the effect of sensing errors on the estimation error. However, the CR 
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bound on the estimation error in the presence of sensing errors cannot be expressed in a simple closed 
form, hence, the estimation error is presented via simulations in Section IVI-CI For analytical tractability, 
uniform sampling is assumed with a constant inter-sample time of T^ seconds. 

We first consider the case where the sensing error can be ignored, i.e., Pf = Pm = 0, thus, Zn = Zn^n. 
The likelihood function of the traffic samples given u assuming perfect sensing is derived in a similar 
manner to ||5] Sec. 6.1] and can be written as 

N-l 
L{z\u) = Y'x{z\u) = Pr(zi|n) TT Y'x{ziJ^i\zi,u) 

i=l 

= u'^ (1 - uf-'^ W Pr,,,,^, {Tc\u), (30) 

where the Markovian property has been applied. Expression (|30l ) can be written as 

xPr^^(r»Pr^f(r,|n), (31) 

where no, ni, n2 and 713 denote the number of (0— s-O), (0— s-l), (1— s-O) and (1— s-l) PU state transitions, 
respectively, from the total of A^ — 1 transitions among N samples. Then, the ML estimator of n, denoted 
by Ura, can be found by solving d\og L{z\u)/du = 0. The value of Um cannot be written in a simple 
closed form and thus, has to be solved numerically. The MSE in Um, denoted by Vu^^n, is lower bounded 
by the CR bound, accordingly, 

K„,JV>Z-'(n,iV), (32) 

where /„ {u, N) is the Fisher information for N collected samples. Specifically, the Fisher information 
is defined as /m(^i -^) = -E' {d log L{z\u)/du) . 

Theorem 2: The lower bound on the MSE in Urn is given as 

(1 - T,) {T, + u- Pen) {FcU -u + 1) 



l;nHn,N) 



Ml + M2 + M-s + M4 + M5 
X u^ (1 - u) , (33) 



where Ml =T'^^XfTc{N-l){l-u)[XfTc{l-u){l+Tc)-2u{l-2u){l-Tc)], M2 = T^^u'^[u{u-l){3N- 
2) + iN- 1)], Ms = -rlu^[u{u - 1){7N - 4) + (27V - 1)], M4 = rcU^[N{5u^ - 5u + 1) + 2m(1 - n)], 
and Afs = Nu^{l - u). 

Proof: See Appendix O D 
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For a fixed observation window, the CR bound for the variance in Um approaches an asymptote as A^ 
increases. This is caused by the correlation between the samples. The lower bound on the CR bound can 
be derived as follows 

lim I^i(n,iV) = !^il^. (34) 



Remarks: Comparing (1341 ) to (1231 ). as A^ approaches infinity, the MSB in estimating u using ML 
estimation with an observation window length of T is the same as that when weighted averaging is used 
with an observation window length of 2T. This is a very useful property for delay-constrained applications 
as the estimation accuracy can be achieved in half the time. The fact that ML estimation requires only 
half the time window length is attributed to the fact that the proposed ML estimator assumes knowledge 
of either Aj or A„, hence, half the information needed to estimate the traffic parameters is assumed to 
be known a priori. The assumption of knowing Aj or A„ beforehand is application specific (the average 
PU off- or on- time may be known to the SU in some network scenarios), where on the other hand, in 
case both Aj and A„ are unknown, joint ML estimation should be used. 

F. ML Estimation of u under Imperfect Sensing 

In the presence of sensing errors, any PU traffic samples vector z can result in an estimated PU traffic 
samples vector z with a non-zero probability. Hence, the likelihood function presented in (|3T| i is modified 
to 

L{z\u) = Y, VT{Zn\u)S{z\Zn), (35) 

n=l 

where Pr(2^„|n) is the probability of occurrence of the PU traffic samples vector Zn and equals the right 
hand side of (|3T| ). and S{z\Zn) is the probability of estimating the PU traffic samples vector as z, when 
the actual PU traffic samples vector equals Zn- S{z\Zn) can be written as 

S{z\Zn) = P™°--(1 - P^)'"!."- 

xP™=— -(l-P^)'"^—-, (36) 

where mo^n,z, 'n^i,n,z, iT^2,n,z, and m3^n,z are the numbers of false alarms, no false alarms, mis-detections, 
and no mis-detections, respectively, that yield the estimated PU traffic samples vector z given the PU 
traffic samples vector Zn. The ML estimator of u can be modified to account for sensing errors and can 
be calculated by solving for the value of u that maximizes the modified likelihood function given in (1351 ). 
The ML estimator as well as the corresponding mean squared estimation error cannot be expressed in a 
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simple closed form and, hence, have to be calculated numerically as shown in Section |Vll 

IV. Estimation of the Primary User Departure Rate A/ and Arrival Rate A„ 

A. Estimation of Xf under Perfect Sensing 

Following 15] Sec. 6.1] the estimation method adopted here is based on ML estimation. Uniform 
sampling is assumed for mathematical tractability. The likelihood function of the traffic samples given 
\f is the same as that in (OTI ) but with replacing the condition on u by a condition on \f. Then, 
the ML estimator of Aj, denoted by Xf, can be found by solving d\og L{z\Xf) / dXj = as in 15] 
Sec. 6.1] yielding A/ = - (V^c) log (- B + ^/W^^IAC\ / {2A) , where A = {u - v?){N - 1), 
B = —2A + A^ — 1 — (1 — n)no — nns, and C = A — uhq — (1 — u)n^. 

B. The CR bound on the MSE in Xf 

The MSE in Xf, denoted by V^ ^, is expressed using the CR bound as in the case of the ML estimation 
of u. Hence, 

V-^^,,>I-'{Xf,N), (37) 



where /„ (A/,iV) = E {d log L{z\Xf)/dXf) 



. Accordingly, we obtain the following theorem. 



Theorem 3: The lower bound on the MSE in Xf is given as 

/-'(A,.iV)= ''<'-r;'F'+""-^'»^-"" . ,38) 

" ' " ' (r,rj^(i-u)(i + r,)(iV-i) 

Proof: See Appendix ID] D 

As for the case of estimating u, the CR bound for the variance in Xf for a fixed observation window 
approaches an asymptote as N increases due to the sample correlation. The lower limit on the CR bound 
can be derived by taking the limit of (|38] | as N approaches infinity, yielding 

C. Estimation of Xn under Perfect Sensing 

The PU arrival rate. An, can be estimated in a similar manner to Xf using ML estimation. The ML 
estimator of A„ can be written as A„ = {I — u) Xf/u where the derivation follows that of A/ and 
is omitted for brevity. It follows that Im{Xn,N) = v?Im{Xf,N)/{\ — u) . Moreover, for a fixed 
observation window, the CR bound for the MSE in A„ approaches A„/ (2Tn) as N increases. 
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Remarks: The expressions for A/ and A„ are functions of u, thus, ML estimation can be used in 
applications where a priori knowledge of u can be assumed. Moreover, the expressions for the MSB in 
\f and A„ are functions of the actual values of Aj and A„, which are not known a priori. Hence, the 
MSB expressions can be used to provide benchmarks for the ML estimation error. Burthermore, the MSB 
expressions can be used to calculate the worst case error in \f and A„, and to show the dependence of 
the estimation error on u, A/, A„, the number of samples, and the observation window length. 

D. ML Estimation of \f and A„ under Imperfect Sensing 

The likelihood function of the traffic samples given A/ or A„ under imperfect sensing can be expressed 
as in (|35] ) with replacing the condition on n by a condition on Aj and A„, respectively. The ML estimators 
for Aj and A„ as well as the corresponding mean squared estimation errors cannot be expressed in simple 
closed forms and consequently, have to be calculated numerically as shown in Section |Vll 

V. Algorithms for the Blind Bstimation of u, A/, and A„ 

In this section, we present algorithms that blindly estimate u, Aj and A„ based on adaptive sampling, 
using the analytical expressions obtained thus far. The assumptions that are necessary for the operation 
of the algorithms are that the off- and on-times of the PU are exponentially distributed. Besides, perfect 
spectrum sensing is assumed, noting that the algorithms can be updated to account for the effect of 
spectrum sensing imperfections. Two algorithms are presented: Algorithm I blindly estimates u assuming 
perfect knowledge of A/ and no a priori knowledge of A„ (see the verbal description in Section IV-AI and 
its summary in Algorithm [!}, and Algorithm II blindly estimates u, Xj and A^ (see the verbal description 
in Section IV-BI and its summary in Algorithm |2]|. 

A. Algorithm \: Blind Estimation of u with Known A/ (or X^) 

Algorithm I is applicable in scenarios where there is a priori knowledge of A/ whereas u and A„ 
are unknown. Note that, the algorithm can be modified to estimate u under the assumption of perfect 
knowledge of A„ with no a priori knowledge of Aj. A practical example for the latter case would be if the 
average on-time of the PU is known (for example, the packet length of the PU follows a certain pattern 
or is fixed) while the rate at which the PU accesses the channel is unknown. The algorithm estimates u 
using the averaging method as described in ^, and the error in Ua is estimated using ^. 

The algorithm operates as follows. The traffic is sampled with an arbitrary initial inter-sample time. 
To, until the sampled traffic state toggles. Then, traffic is sampled for an arbitrary initial number of 
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samples, A'o, with the inter-sample time Tq. Note than increasing Tq and Nq increases the accuracy 
of the initial estimate, but, on the other hand, increases the estimation delay. For the remainder of the 
algorithm, the inter-sample time is determined while taking the correlation between the traffic samples 
into consideration. Corollary |2] provides an expression for the expected decrease in the MSE in Ua for 
each additional sample. The maximum decrease in the MSE in Ua is given by 

Tn—^oo 

^ K„,;v(2iV + l)-n(l-n) 

{N + 1)2 • ^ ""^ 

Accordingly, the parameter D^^^jy^i can be expressed as D^^ ^r+i = -D™^^^^ — DqTn, where Dq = 
(■jv+ip Si=o Ilk=N-i ^fc- Parameter -D^^^^^+i is independent of T/v, thus, varying T/v can only affect 
the DqTn term in Z)^^ ^v+i- Setting T/v to ^ is equivalent to multiplying Dq by a factor of e~". The 
parameter a is chosen so that a compromise between the reduction in the estimation error per sample 
and the delay in taking the new sample is reached. For lower values of a, the expected decrease in the 
MSE in Ua per sample is lower, while at the same time, the average delay in taking the new sample is 
lower. 

There are three different conditions for terminating the algorithm that are used depending on the 
application. The algorithm may terminate when (i) a predetermined number of samples are taken (energy- 
constrained applications), or (ii) after a certain observation window length is reached (delay-constrained 
applications), or (iii) when a target expected estimation error is reached. For an energy-constrained 
application where the sensing energy budget, and hence the total number of samples that are to be taken, 
is limited, the first termination condition is used where the total number of samples equal Nth- On the 
other hand, for a delay-constrained application where the total observation window is bounded by a time 
threshold, T^^, the second termination condition is applied. Finally, the third termination condition is 
used if the application involves estimating u with a target average estimation error Vth- To ensure that 
the target average error is always reached, the expected error is calculated using dS) for u G [0, 1], using 
the operating values of N, Aj, and T. The algorithm is terminated if the worst case expected estimation 
error, calculated over all possible values of u, is less than Vth- Note that the value of u that yields the 
highest estimation error cannot be known a priori as described in Section [VI-AI The proposed algorithm 
is summarized in Algorithm [T] 
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Algorithm 1: Algorithm I: Blind Estimation of u with known A/ (or A„) 

Require: Tq, Nq, a 

Ensure: T < Tth or N < Nth or max{T4^^7v} < V^/j (user specified) 

1: while Ua = or Uq = 1 do 

2: take a next sample after Tq s 

3: calculate Ua using ^ 

4: if T > Tth or N > Nth (if target Tth or Nth is required) then 

5: Terminate 

6: end if 

7: end while 

8: while iV < A'^o do 

9: take a next sample after Tq s 

10: calculate Ua using (|2ll 

11: if T > Ti/j or A^ > Nth (if target Ti/^ or Nth is required) then 

12: Terminate 

13: end if 

14: end while 

15: while T < Tth or N < Nth or max{T4^ at} > T4/, (user specified) do 

16: calculate Tjsi = ^a 

17: Take a next sample after T/y s 

18: calculate Ua using Q 

19: calculate T4„,Ar for all u G [0, 1] using (|5]) (if target Vth is required) 

20: end while 



B. Algorithm II; Blind Estimation of u, A/ ancf An 

Algorithm II is directed for scenarios where there is no a priori knowledge of u, Xj, and A„. The duty 
cycle is estimated using the averaging method as described in ([2]). Parameters Aj and A„ are estimated 
using ML estimation as described in Section |IVl hence, unlike Algorithm I, the inter-sample time is 
kept constant at Tq. The algorithm terminates when a target average estimation error in u, denoted by 
V„ j/j, and a target average estimation error in A/- or A„, denoted by Vx^th^ are reached. To ensure that 
the targeted estimation errors are met, the estimation errors are calculated at all algorithm iterations for 
the full range of u, and Aj or A„. Moreover, it follows from Section |IV] that the asymptotic lower bound 
on the eiTor in Xf is greater than that in A„ for n > ^ and approaches infinity as u tends to 1, while the 
asymptotic lower bound on the error in A„ is greater than that in Xf for n < i and approaches infinity 
as u tends to 0. Thus, the algorithm is designed to terminate when Ua and A/ reach the target average 
estimation error if Ua < ^, and when Ua and A„ reach the target average estimation error if Ua > ^■ 
Besides, The lower bound on the estimation error for Xj and A^ is proportional to Xf and A„ as shown in 
Section |IVl Accordingly, the error in A/ and A^ cannot be guaranteed to meet specific targets unless A/ 
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Algorithm 2: Algorithm II: Blind Estimation of u, A/, and A^^ 



Require: Tq 

Ensure: max{T4„,7v} < K,ih and maxjV^^^^} < Vx,th for Ua < 0.5, and max{T4^,Ar} < Vu,th and 
max{ys^^ ^} (MSE in A„) < Vx,th for «« > 0.5 



1 

2: 
3: 
4: 
5: 

6: 
7: 
8: 
9: 
10 
11 



while Uq = or {ta = 1 do 

take a next sample after Tq s 

calculate ita using ^ 
end while 

while {u<j < 0.5, and max{T4^^Ar} > K,ih or maxjV^^^^^} > Vx,th} or {ua > 0.5, and 
max{T4^,Ar} > Vu,th or maxj^s^^^^} > Vx,th} do 

take a next sample after Tq s 

calculate Ua using Q 

calculate Vu^^n for all n G [0, 1] and for Xf = Amin using (|5]l 

calculate A/ and A^ following Section ITVl using Ua, N, uq, and n^ 

calculate V^ ^ and V^ ^ for all u G [0, 1] and for A/ = A„ = Amax following Section |IV] 
end while 



and A„ are upper bounded. Moreover, we assume that A/ and A„ are greater than zero, as otherwise, the 
PU would be either always on or always off. Thus, in this section we assume that A/, A„ G [Amin, Amax]- 
Furthermore, since the error in estimating u increases with decreasing Aj as presented in (JS), the worst 
case estimation error in u is calculated while substituting Xf = Amin- On the other hand, as the estimation 
error in Xj and A„ increases monotonically with Xf and A„, respectively, the worst case estimation error 
in Xf and A^ is calculated while substituting Xf = A„ = Amax- The proposed algorithm is summarized 
in Algorithm |2] 

VI. Numerical Results 

We present numerical results, confirming the theory provided in Section |llll that compare the accuracy 
of the different estimation methods of u. Moreover, we show the accuracy and the asymptotic behavior 
of the ML estimation of Xf, based on the theory obtained in Section JV] Furthermore, we present results 
for the estimation error in u and Xf under spectrum sensing imperfections. Note that the performance 
of the estimation of A„ is similar to that of Xf, and thus is omitted to eliminate redundancy. Next, 
we present results showing the performance of the proposed blind estimation algorithms. Besides, we 
validate the correctness of the developed mathematical expressions through comparison with simulation 
results developed in Matlab version 7.10.0.499. Note that, the root mean squared (RMS) error is used as 
a metric for quantifying the estimation accuracy instead of the MSE for convenience. The RMS error is 
simply calculated as the square root of the MSE. 
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In the numerical evaluation, as an example, typical values for the PU traffic parameters are used 
following the results in Q, lITSl . which are representative sources of information of the temporal 
utilization of radio resources in popular radio access systems. In specific, we focus on GSM 1800 
downlink traffic, as it is the only service common to both studies, with detailed parameters found in lITS] 
Tab. VIII] and 121 Tab. 3]. We therefore assume that the PU duty cycle, u, is in the range of [0.3,0.6] 
(0.30 in HH and 0.62 in ||2l), whereas the PU departure rate. A/, is in the ranged of [0.4,0.9] s~^. 

A. The RMS Error in Estimating u 

1) The Variation of the RMS Error with the Number of Samples: The relationships between the total 
number of samples, N, and the RMS error in the estimate of u are plotted in Fig. |2]for (i) the averaging 
estimator with non-uniform sampling using the optimal inter-sample time sequence, (ii) the averaging 
estimator with uniform sampling, (iii) the weighted averaging estimator using the optimal weighting 
sequence, and (iv) the ML estimator The RMS error in estimating u, denoted by Ru,n, is plotted for 
an observation window of T = 50 s where N is increased from 40 to 150 samples, which represents a 
typical size of sample sets in traffic sampling. The variation of Ru^n with N ior u = 0.3 and u = 0.6 is 
shown in Fig. |2(a)| and Fig. |2(b)[ respectively. The PU departure rate, Xj is set to 0.4 s~^ and 0.9 s~^. 

The results show that ML estimation outperforms all averaging based estimation techniques where the 
resulting RMS error can be reduced by up to 24%. This is because the proposed ML estimator assumes 
a priori knowledge of A/-. Moreover, the optimized averaging estimator with non-uniform sampling and 
the optimized weighted averaging estimator yield almost the same estimation error, with a narrow margin 
below the averaging estimator with uniform sampling. Besides, for the same u, higher A/ yields lower 
estimation error due to the reduced sample correlation. Furthermore, the figures emphasize the fact that the 
estimation error reaches an asymptotic value as N is increased. Finally, all results are verified via Matlab 
simulations where the simulation results match the theoretical expressions except for ML estimation 
where the simulation-based error is higher than the theoretical expression. This is because the CR bound 
provides a lower bound on the error that is attained asymptotically as N increases. 

2) The Asymptotic RMS Error: The RMS error in the estimate of u as A^ tends to infinity reaches an 
asymptote as shown in Section |llll Fig. [3] shows the relationship between the asymptotic RMS error in 
the estimate of u, denoted by Ru^oo, and the observation window length, T, for the different estimation 

'As the measurements in (2] were performed in the discrete domain, we directly converted the parameters of all geometric 
distributions given in |J2) Tab. 3] into exponential distribution parameters, where the minimum and maximum reported PU 
departure rate is 0.48 s^^ and 0.9 s^^, respectively. Note that the arrival and departure rates are not reported explicitly in tl8i . 
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Fig. 2. RMS en'or in tiie estimate of u as a function of tlie number of samples A'^ for T = 50 s and different values of it. The 
RMS en^or is plotted for the four estimation methods; US: Averaging with uniform sampling, NS: Averaging with non-uniform 
sampling, WS: Weighted averaging, and MLE: Maximum likelihood estimation. Simulation results (Sim.) are plotted to verify 
the mathematical model (An.). 
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Fig. 3. Asymptotic RMS error in the estimate of u, Rn.oa, as a function of the observation window length, T. Ra.ao is 
plotted for the four estimation methods; US: Averaging with uniform sampling, NS: Averaging with non-uniform sampling, WS: 
Weighted averaging, and MLE: Maximum likelihood estimation. The traffic parameters used: {u = 0.3, A/ = 0.9 s^^}, and 
{u = 0.6,Xf =0.4s-^}. 



techniques. Two different sets of PU traffic parameters are used; u = 0.3 and Xf = 0.9 s^^, and u = 0.6 
and Xf = 0.4 s~^. Ru^oo for the averaging estimator with non-uniform sampling is calculated numerically 
since a closed form expression is not available. The results show that Ru^oo for ML estimation is lower than 
that for the other estimation techniques, and the performance of the weighted averaging estimator and the 
averaging estimator with non-uniform sampling is almost identical and surpasses that of the averaging 
estimator with uniform sampling. Moreover, as proven in Section IIII-E[ Ru^oo for ML estimation is 
identical to that for weighted averaging estimation if the observation window length, T, is doubled. 

3) The Variation of the RMS Error with the Duty Cycle: The relationship between Ru^n and u, is 
presented in Fig. ID The plot compares the error for the averaging estimator with uniform sampling with 
that of the ML estimator. For this setup, u is increased from to 1 while T and N are kept constant at 100 
seconds and 100 samples, respectively. The error is presented for different values of Xf representing traffic 
with different levels of correlation. ML estimation achieves a more accurate estimate than averaging-based 
estimation, yet the gap in performance decreases with higher A/. The value of u which results in the 
highest estimation error is greater than 0.5 and approaches u = 0.5 with increasing Xf. The skew in the 
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Fig. 4. RMS error in the estimate of it as a function of u for TV = 100 and T = 100 s. The RMS error is plotted for averaging 
with uniform sampling (US) and the CR bound is plotted for maximum likelihood estimation (MLE). The traffic parameters 
used: A/ G {0.1, 0.5, 1, 2} s~^. Simulation results (Sim.) are plotted to verify the mathematical model (An.). 



figure is attributed to the error added by sample correlation, which increases with u for the same Aj. For 
higher A/, the effect of sample correlation decreases and the mean squared estimation error approaches 
j^^ , as explained in Section ITlI-Al[ which is symmetric in u and is maximized at u = 0.5. Again, the 
results are verified via simulations where the simulation results match the theoretical expressions except 
for ML estimation as the CR bound provides a theoretical lower bound on the error. 

Key Message: ML estimation is recommended for estimating u when a priori knowledge of Xf 
or A„ is available. Moreover, optimized averaging under non-uniform sampling and optimized weighted 
averaging yield almost the same estimation error, and result in a lower estimation error than averaging 
under uniform sampling. 



B. The RMS Error in Aj under Uniform Sampling 

The PU departure rate, Aj, is estimated using ML estimation and hence, the estimate accuracy is lower 



bounded by the CR bound. Fig. |5(a)| shows the square root of the CR bound as well as the RMS error in 
Xf, obtained by simulations, for different traffic parameters and T = 50 s, as a function of the number 
of samples, N. The estimation error reaches an asymptote as N is increased for a fixed observation 
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Fig. 5. CR bound [5(a)l and asymptotic CR bound [5(b)l on the error in A/. The used traffic parameters: {it = 0.3, A/ = 0.4 s~^}, 
{u = 0.3, A/ = 0.9 s~^}, {u = 0.6, A/ = 0.4 s }, and {u = 0.6, A/ = 0.9 s~^}. Simulation results (Sim.) are plotted to 
verify the mathematical model (An.). 



window length, which, again, is attributed to sample correlation. The asymptotic value of the CR bound 
was derived as ( [39l ) and its square root is presented in Fig. |5(b)| where it is clear that the asymptote 
decreases with increasing the observation window length. The analytical result in (|39l ) is verified through 
simulations where N was set to a value that is higher than TXf/u by orders of magnitude. 

Key Message: The PU departure rate can be estimated using ML estimation where the estimation 
error is lower bounded due to sample correlation. The lower bound on the error can only be decreased 
by increasing the total observation window length. 



C. RMS Error in Estimating u and A/ under Sensing Imperfections 

The impact of sensing imperfections on the estimation of u and Aj is presented in Fig. |6(a)| and 
Fig. |6(b)[ respectively. The assumed parameters are A^ = 0.9 s~^, u = 0.3, and T = 50 s. The estimation 
error is shown for Pf = Pm = 0, Pf = Pm = 0.05, and Pf = P^ = 0.1. The estimation error under 
sensing imperfections for the averaging estimator of u is expressed in closed form in (|T9] l. and the 
error for the weighted averaging estimator using the optimal weighting sequence is expressed in closed 
form in (|29l ). Analysis and simulation-based results are plotted for both estimators (the special case of 
uniform sampling is used for the averaging estimator), where the simulation results match the theoretical 
expressions for both estimators. On the other hand, the impact of sensing errors on the error in the ML 
estimators of u and A/ is only expressed via simulations using the modified likelihood function expressed 
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in (l35]) . 

Key Message: For estimating u, the ML estimator outperforms the averaging estimators under sensing 
imperfections. The effect of sensing errors on the error in estimating A/^ is noticeable where the RMS of 
the estimation error increases by up to 91% for Pf = Pm = 0.1, compared to perfect spectrum sensing. 
Finally, the impact of the sensing imperfections on the PU traffic parameters estimation error decreases 
as N increases. 

D. Algorithm I: Performance of the Proposed Duty Cycle Estimation Algorithm for Constrained N 

This section presents the performance of Algorithm I when the total number of samples is constrained. 
A typical application would be traffic estimation by an energy-constrained node. The initial number of 
samples, A'^o. is set to 5 samples. The initial inter-sample time, Tq, is set to two different values, 1 and 
10 seconds, to investigate the performance of the algorithm under different initial conditions. After the 
initial Nq samples, the inter-sample time is adapted as described in Section |Vl Three different values 
of a G {1,2,5}, are selected to show the compromise between the estimation accuracy and the total 
observation window length. The performance of the algorithm is compared to that of uniform sampling 
with Tu set equal to the 2 different values of Tq, i.e., the algorithm is compared to the case where the 
inter-sample time is kept constant at the initial conditions without adaptation. The traffic parameters are 
u = 0.6 and A/ = 0.9 s~^. The RMS error in the estimate of u, Ru^n, is presented in Fig. 7(a) and the 



equivalent total observation window length is presented in Fig. |7(b)| As a reference, the theoretical lower 
bound for the averaging-based estimation error as T tends to infinity is also plotted. 

The results show that the algorithm can blindly achieve an RMS estimation error that is only 1.5% 
higher than the theoretical lower bound for q = 5, Tq = 10 s, and Nth = 100 samples. Uniform sampling 
alone can achieve a low estimation error but at the expense of a notable increase in the total observation 
window length. Moreover, the algorithm has lower dependence on the initial conditions compared to 
uniform sampling, as the algorithm adapts the inter-sample time according to the traffic parameters. 
Furthermore, higher a yields a smaller estimation error but causes an increase in the estimation duration, 
thus, a can be tuned according to the application constraints. 

Key Message: The proposed algorithm can blindly estimate u for a constrained N by adapting the 
inter-sample time according to the estimated u. The compromise between the estimation eixor and the 
total observation window length can be controlled by tuning the algorithm parameter a. 
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and T = 50 s. The RMS error is plotted for three estimation methods; US: Averaging with uniform sampling, WS: Weighted 
averaging, and MLE: Maximum likelihood estimation. Simulation results (Sim.) are plotted to verify the mathematical model 
(An.) where applicable. 
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Fig. 7. The performance of Algorithm I for constrained Nth, Fig. |7(a)| The RMS estimation error in u as a function of 
the number of samples threshold Nth ', Fig. |7(b)| The total observation window length as a function of the number of samples 
threshold Nth- The estimation error and total observation window length for uniform sampling are plotted for comparison. The 
theoretical lower bound on error is plotted as a reference; Alg.: Algorithm I, Uni.: uniform sampling. The used parameters: 

u = 0.6, \f = 0.9s-\ iVo = 5, To G {l,10}s. 
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Fig. 8. The performance of Algorithm I for constrained RMS error. Fig. |8(a)| The achieved RMS estimation error in u as 
a function of it; Fig. |8(b)| The average total number of samples as a function of u; Fig. |8(c)| The average total observation 
window length as a function of u. The used parameters: A/ — 0.9 s^^, No — 50, To = 50 ms. 



E. Algorithm I: Performance of the Proposed Duty Cycle Estimation Algorithm given a Target Estimation 
Error 

In this section, the duty cycle is blindly estimated until the RMS error reaches a targeted value using 
Algorithm I. The target RMS estimation error is 0. 1 and the algorithm is tested for u ranging from 0. 1 
to 0.9, while \f is set to 0.9 s~^. The algorithm parameters are set to Nq = 50 samples, Tq = 50 ms, 
and a G {1, 2, 5}. The RMS estimation error, the average total number of samples, and the average total 
observation window length are plotted in Fig. 8(a) Fig. |8(b)| and Fig. 8(c)[ respectively. 

Key Message: The achieved RMS estimation error is always reached for all values of u. The reached 
error is less than the targeted error for most cases as the algorithm targets the worst case error Furthermore, 
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the results emphasize the compromise between the number of samples and the observation window length 
where a can be chosen according to the sensing energy and sensing delay constraints. 



F. Algorithm II; Joint Estimation of u, and Xj (or A„j for a Target Estimation Error 

The targeted RMS estimation error in u is set to 0.1, and the targeted RMS error in Xj or An is set to 
0. ls~^. The algorithm is tested for u ranging from 0.1 to 0.9, and A/ G {0.1,0.5,0.9} s~^. Parameters 
Amin and Amax are set to 0.1 s^^ and 1 s~^, respectively, and the inter-sample time, Tq, is set to 50 ms. The 
RMS estimation error in u is presented in Fig. |9(a)[ and the RMS error in A/ or An is shown in Fig. |9(b)| 
The results show that the constraints on the estimation error are reached for all of the tested values of u 
and A J. The estimation error in u is higher for lower Aj due to the increased sample correlation, while 
the error in A/- increases with Xj which complies with Section ITV-B I Furthermore, since the inter-sample 
time is constant and the algorithm terminates when the worst case error (covering the full considered 
range of the traffic parameters) matches the target error, the total observation window length and number 
of samples are constant for all considered values of u and Xf. For this specific setup, the total observation 



window equals 290 seconds, which is higher than that reported for Algorithm I, as shown in Fig. |8(c) 
where Xf was assumed to be known a priori. 

Key Message: The proposed algorithm can blindly estimate all traffic parameters while satisfying a 
targeted estimation error. 



37 

hn,i,j = Z^a=l llfc=l ^k, i = j = ^i 

"-rijij — 1 rLn,i,j ^ llfc=l-'-fc) ^ ^ J-jJ > i, 

h'n,i,j = ■^ '^n,i,j + "■n,i-l,j — 1 ~ h,n,i-l,j ~ h,n,i,j — l = ~ llfc=i 's' J > "* > Ij 

/in,i,j + hn,i-l,j-l - hn,i-l,j - K^iJ-l = Y.a=i \\.k=i ^fe + YJ'a=l \l\=a ^fc, J = « > 1, 

/inji, Otherwise. 

(41) 



VII. Conclusions 

In this paper, we developed a mathematical framework that quantifies the estimation accuracy of PU 
traffic parameters in the form of the mean squared estimation error. We derived the Cramer-Rao bound on 
the accuracy of estimating the mean PU duty cycle, u, using uniform sampling under perfect knowledge of 
either the mean PU off- or on-time. We also analyzed the estimation error in u for a variety of estimators 
based on sample stream averaging and maximum likelihood estimation. We proved that for the equal- 
weighted averaging-based estimation of u, the estimation error is convex with respect to the inter-sample 
time sequence and we derived the optimal sequence and the corresponding error. Moreover, we derived 
the optimal weighting sequence when weighted averaging is employed. Furthermore, we showed that 
the maximum likelihood estimation of u outperforms all averaging based estimators, provided that a 
priori knowledge of the mean PU off-time or the mean PU on-time is available. Regarding the mean 
PU off- and on-times, we formulated the estimation enor bounds, for all unbiased estimators, in the 
form of the Cramer-Rao bounds. We showed that the estimation error for all PU traffic parameters for 
a fixed observation length is lower bounded due to sample correlation, where the bound can be reduced 
by increasing the total observation window length. Besides, we demonstrated the impact of spectrum 
sensing errors on the estimation accuracy for all PU traffic parameters, including analytical results where 
applicable, and showed that the effect of sensing errors on the estimation accuracy of u can be eliminated 
by increasing the number of traffic samples. Finally, we proposed algorithms, based on the derived error 
expressions, for the blind estimation of u for constrained number of samples, observation window length, 
and expected estimation error. We concluded the paper by proposing an algorithm for the joint estimation 
of all traffic parameters that successfully achieves a target estimation error 
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Appendix A 
Proof of Vu^,,^n 
Proof: The expression presented in (|T9] l can be proved by mathematical induction. First assume, 
without loss of generality, that I — Pj — P^ 7^ 0. Define P^ = Pr(i: = Zi\T), i.e., the probability of 
observing the estimated PU traffic samples vector Zi given T. Moreover, define P as a vector containing 
all values of i-*„ such that P„ is the nth element of P. Furthermore, define 5 as a vector where Si, 
defined in Section Ull-Bi is the fth element of S. Hence, Vu^ ^^n = X]i=i ^f-^i ~ '^■ 

For the base case with N = 2, Z = [00,01,10,11] and S = -1/(1 - Pf - Pm)[-Pf,-Pf + 

l/2,-Pf + l/2,-Pf + l].Besides,P=[{l-u)Pf{ProoPf + PmPm)+uPmiFnoPf + PniPm),i^- 
u)PfiPm Pf + Proi P,n) + nP„(Prio Pf + Prn P„), (1 - n)P^(Proo Pf + Proi Pm) + nP,„(Prio Pf + 
Prn Pm), (1 - n)P/(Proo Pf + Proi An) + uPm(Prio Pf + Pm Pm)], where Pf = 1 - Pf and P^ = 
1 — Pm- Accordingly, 

u{l - u)Tc , u(l - u) 

Vu„ ..2 = ^ \- 



2 2 

, nPn.(l-P>n) + (l-^)P/(l-P/) 
+ 2(1 - P/ - P.)2 ' ^"^-^^^ 

which equals ( fT9l ) for N = 2, hence, proves the base case. 

Showing that (fT9l) holds for A^ + 1 while assuming that it is true for A^ is sufficient for proving (IT9] l. 
Subscripts A^ and A^ + 1 are added to P„ and Sn to differentiate between cases with A^ and A^ + 1 
samples. For A^ + 1 samples, V^^ ,,N+i = z2i=i 'SiN+iPi,N+i — u and can be expressed as 

N^Vu N (2N + l)n2 



where 



2W+1 / jv \ 



i=l \n=l 



and 



2 

^ (AT + 1)2(1 -P^-P„,)2 (^-44) 



2N + 1 



f^2 = 2^ Pi,N+l{Zi^]\f+i — Pf) 



i=l 



^ (AT + 1)2(1 -P^-P„)2- (^-45) 
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Note that ^j tv+i in ^i and ^^2 denotes the estimated traffic sample number A^ + 1 in the traffic sample 
vector Zi. Oi is a recursive expression that can be simplified to 

where the proof is omitted for brevity. Moreover, ^^2 can be simplified to 

Finally, substituting ( fT9l ) for Vu^ ^^n in (IA.43I) and using the simplified expressions for Q.\ and Vl^, we 
obtain 

2m(1 - u) ^-^ V^ -TT -^fc^/ 'u(l-'u) 

uP^{l-P^) + {l-u)Pf{l-Pf) 
{N+l){l-Pf-P^)^ 

This corresponds to (|T9] l with the number of samples set to A^ + 1, hence, proves (|T9] l by mathematical 
induction. D 



Appendix B 
Proof of Convexity of Vu^^n {T) 
Proof: In this section, expression dS) is proved to be convex by showing that the Hessian of Vu^^n (T), 
V^T4„,Ar (T), is positive-semidefinite OTI . Denote the Hessian of Vu^^n (T) hy H = [/ia,fe]- -ff is a sym- 
metric {N - l}-by-{7V - 1} matrix and can be expressed as ha,b = ^"^^^7"^ i^j Ya=i Y!j=b I^k=i ^fc' 
Va < h. Since ' )y7" (if ) ^0' showing that H is positive-semidefinite proves that H is positive- 
semidefinite, where if = 2u(i-u) ( X" ) H. H is, proved to be positive-semidefinite by showing that all 
of its A^ — 1 leading principal minors are non-negative. Define Hn = [hn,i,j\ as the upper left n-by-n 
comer of H. Showing that the determinant of if„ is non-negative for n E {1, 2, • • • A^ — 1} proves that 
H is positive-semidefinite. Define Hn = \hn,i,j\ as a symmetric matrix where the elements of Hn are 
defined in (|4T]) . given at the top of the previous page. 

It is easy to show that | Hn \ = \ Hn \ as H^ is formed by performing row and column addition operations 
on Hn and an even number of sign changes. Moreover, Hn is diagonally dominant, that is, for every 
row of the matrix, the magnitude of the diagonal element is greater than or equal to the summation of 
the magnitude of the other non-diagonal elements: \hn,i,i\ — Ylj^i \hn,i,j\ = Y,a=n YYk=i^k > 0,Vi G 
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{!,••• ,n},Vn S {I,-'' ,N — 1}- A real symmetric diagonally dominant matrix with non-negative 
diagonal elements is positive-semidefinite BOl Ch. 6]. Hence, \Hn\ > 0, and accordingly, \Hn\ > 0, for 
n G {1, 2, • • • A^ — 1}. Consequently, H is positive-semidefinite and Vu^^n (T) is convex. D 

Appendix C 
Proof of the Lower Bound on 14„,Af 
Proof: We first start by simplifying the expression for the Fisher information. Let ^n = d log L{z\u)/du 
where A^ samples are used for estimation. Note that Im {u, N) = E [^%] ■ Using ([T]) and (OTT ). <l>jv can 
be written as 



Zl- U ^0 



U{1 — u) 



+ $1 



ni no 



Proi(rc|n) Proo(rc|u) 



(C.49) 



Pni{Tc\u) Prio(rc|u) 

where $o = (1 - r^) - XfTcTc/u, $i = (1 - r^) + X/TcTc (1 - u) jv?, and P^ = e'^i^'l'^. We now 

apply mathematical induction to prove (1331 ). Starting with the base case of A^ = 2, /„ (u, 2) = E [^2] . 

The expectation is evaluated over all four possible sample sequences. Thus, Im. (u, 2) = X]j=i 4>2 i Pr(z = 

Zi\T), where (j)j\f^i is ^n evaluated for Zi. Furthermore, Z = [00,01,10,11] and the corresponding 

values of Pr(z = Z,\T) = [(1 - u)Fm{Ti), (1 - n) Proi(Ti), nPrio(ri), txPrn(ri)], where Fr,y{-) 

is as defined in ([T]). Denote the vector of all values of cl)2,i for i G {1,2,3,4} by (/)2- Hence, 02 = 
1 _ u{i-r^)-\fT,r, 1 u{i-r,)-\fT^r, i _ u'^ (i-r,)+\fT,r^{i-u) 1 u^{i-r,)+XfT,T,{i-u) 

(m-I) MProo{Tc|M) '(«-!)"'" «Proi(Te|M) ' m «2Prio(Te|u) ' m + u'^ Prii{Ta\u) 

Accordingly, Im {u, 2) can be written as 

M2,l + M2,2 + M2,3 



-fm (u, 2) 



(1 - Pc)(rc + U- ure)(uPe - "U + 1) 
1 



(C.50) 



M'^(l — u) ' 

where M2,i = A/rcr2(l-n) [A/rc(l - 'u)(l + Pc) - 2u(l - 2n)(l - P^)], M2,2 = P^n^ [iu{u - 1) + 1]- 
r^u^ [10u{u - 1) + 3], and M2,3 = PcU^ [2{5u^ - 5u + I) + 2u{l - n)] + 2n^(l-n). Expression (IC30l) 
is equivalent to (1331 ) for N = 2, which proves the base case. Assuming that (|33] | is true for any N, proving 
that (|33] ) holds for A^ + 1 completes the proof. The Fisher information for A^ + 1, can be expressed as 

2N + 1 



i=l 



X Pr(z(^+i) = zf +^)|r(^+i)), (C.51) 
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where the superscripts (N) and {N + 1) are appended to z, Zn, Z, and T, to differentiate between the 
cases with A'^ and A^ + 1 samples, respectively. The set of all possible sample sequences of length N, 
Z^^\ can be split to two subsets, ^(^'°) and Z^^'^\ which represent the set of all sample sequences 
ending with and 1, respectively. Thus, the summation over the set z(^+^) presented in (IC.Sll l can be 
split to summations over the sets Z^^'^^ and Z^^'^^ as follows 



I„^{u,N+l) 



E P^^.^ [ [^^^^ - Prootr,|n) ) P^oo(7;|n) 



2<"'g^(Af,0) 



+ (t>NA + 



$0 



Pm{Tc\u) 



Proi(Tch) 






2}"'g2(".i) 



+ (pNA + 



$1 



Prii(Teltx) 



Prii(re|7x) 



E 



Pri,jv ^' 



02 



2UV)g^(JV,0) 



u{l-Tc){uT,-u + l) 



^ . E (l-r,)(l 



Pr, TV ^ 



l2 



2}"'g2(iv,i) 



(i-re)(i-u)[u + (i-n)r. 



+ -^m(w,iV) 



(^r, -u + i)[ii- TcW + A/rje(i - u)] ^ 
#(1 -«)(!- rc)(rc + u- nrc)(urc - -u + 1) 



+ 



(1 - u)^Tc + u- uT^) [u(l - Pe) - A/T,r, 



f-^a-ci 



u^{l - u){l - Tc)iTc + u- uTc)iurc -u + l) 

+ Im{u,N), 



(C.52) 



where Pr^.y = Pr(2^ = Zx\T^). The expression for Lin{u,N + 1) given in ( IC.52I ) can be shown to 
match that in (1331 ) for A^ + 1 by substituting Im{u,N) in (IC.52I ) by (1331 ) and simplifying the resulting 
expression. This concludes the proof. D 



42 

Appendix D 
Proof of the Lower Bound on V^ ^ 

Proof: As in Appendix O let Gat = (91ogL(z|A/)/(9A/, where Im (A/,iV) = E [6^]. Using ©, 
and (|3T| ) (with the condition on u replaced by a condition on Aj), ©tv can be written as 



e 



N 



U 

+ (1 



ni no 

"" Proi(r,|A/) Proo(T,|A/) 

"-2 ^3 



(D.53) 



^Prio(re|A/) Prn(r,|A/) 
We now apply mathematical induction to prove (|38] ). Starting with the base case of A^ = 2, 1^ (A/, 2) = 
E [©1] — Si=i ^2 j -P'^(-^ ~ ^i\T), where 0^,1 is ©at evaluated for Zi, and Z and the corresponding 
values of Pr(z = Zi\T) are as defined in Section [Cl Denote the vector of all values of 6*2,1 by 62, hence, 

''^ - [Proo(T,|A^)' Proi(r,|A;)' «Prio(r,|Aj)' «Prn(T,|A^)J • ^ccormngiy, 

'- ^^^' '^ - n(i-r.)[r. + .(i-r.)2(i-.)]' (°-^4) 

which is equivalent to (|38] ) for N = 2. Assuming that ( [38] ) is true for any N, proving that (|38T l holds for 
A^ + 1 is sufficient for proving that (|38] | is valid for any A^. For A^ + 1, the Fisher information can be 
expressed as 

2JV+1 

X Pr(z(^+i) = zf +^)|r(^+i)), (D.55) 

where the superscripts (A^) and (A^ + 1) are appended to z, Zn, Z, and T, to differentiate between the 
cases with N and A^ + 1 samples, respectively. As in Appendix O the summation in (ID.55I ) can be split 
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to summations over sample sequences of length N ending with and 1 . Thus 



I„(A/,7V + 1) 



2(N)g2(Af,0) 



ON,i - i„„(yjr.) ) Pm{Tc\u 



+ 9m i + 



TcTc 



u{l-T,) 



+ 5Z P^^.^ 

2(™)g2(JV,l) 

+ Ona + 



u^ + u(i - u)r. 



Proi(rc|n) 



^N,i + ^:^r^) PnoiTc\u 



Prii(T,|u) 



Substituting ([T]) in ( ID.56I ) and simplifying the resulting expression 

+ \" Pr- ^^ (rcT.)^(l-«) 

+ Z^ ^^hN u^i-r,)(r,+u-r,u) 
2(™)e^(jv,i) 

+ /m(A/,iV) 



n(i-re)[re + «(i-re)2(i-^)] 

+ /^(A/,iV). 



Finally, substituting dM) in (ID.57I) yields 

I;„(A/,iV+l) = 



iv(r,T,)2(i + re) 



(i-r,)[re + tx(i-r,)2(i-n)] 

l-n 



n 



This proves that (1381 ) holds for A^ + 1 and thus concludes the proof. 



(D.56) 



(D.57) 



(D.58) 



D 
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